Effect of discontinuities in Kohn-Sham-based chemical reactivity theory 
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We provide a new derivation of a formula for the Fukui function of density-functional chemical 
reactivity theory which incorporates the discontinuities in the Kohn-Sham reference system. Orbital 
relaxations are described in terms of the exchange-correlation (XC) kernel, i.e., the derivative of 
the XC potential with respect to the density and it is shown that in order to correctly measure 
the reactivity toward a nucleophilic reagent a discontinuity of the XC kernel has to be taken into 
account. The importance of this finding is illustrated in model molecular systems. 



I. INTRODUCTION 



Predicting how molecules respond to external pertur- 
bations is an important subject in theoretical chemistry. 
On a fundamental level this entails a very difficult prob- 
lem as molecules are composed of interacting electrons 
and nuclei which require a solution to the many-body 
Schrodinger equation. In an attempt to simplify the 
description a set of reactivity descriptors and associate 
empirical equalization principles have been formulated, 
constituting what is known as chemical reactivity the- 
ory (CRT).^"^"^ An approach which naturally combines 
with CRT is density functional theory (DFT),"'i^ which 
is an exact framework for treating the electron-electron 
interaction in terms of only the electronic density. Nor- 
mally, the Kohn-Sham (KS) formulation^^ is used in 
which the density is calculated from a fictitious system 
of non-interacting electrons moving in an effective local 
multiplicative potential - the KS potential. It is well 
known that an independent-particle description using a 
local potential sometimes introduces singularities in the 
KS quantities and that many physical properties are cru- 
cially dependent on this singular behavior. A classic ex- 
ample is the dissociation of closed-shell molecules com- 
posed of open-shell atoms. In order to dissociate with 
a correct integer number of electrons on each atom a dis- 
continuous positive step in the exchange-correlation (XC) 
part of the KS potential has to develop over the atom 
with the larger ionization potential. Important progress 
was made when it was understood that this feature can 
be related to a derivative discontinuity of the total energy 
as a function of particle number. In order to make this 
identification DFT was generalized to ensembles allowing 
for fractional charges. 

In CRT the use of non-integer number of particles is 
at the basis of all definitions of reactivity indices since 
the reaction probability is mostly measured in terms of a 
sensitivity to a global or local change of particle number. 
In KS-based CRT one therefore expects singularities as 
those discussed above to become important. The aim of 
this paper is to show an example which points out this 
fact, namely in the determination of the Fukui function, 
a local reactivity index. In this case, the quantity having 
the crucial discontinuity is the XC kernel defined as the 
functional derivative of the XC potential with respect to 



the density. The XC kernel is the quantity that accounts 
for the effects of orbital relaxation, which may produce a 
large difference when predicting the reactivity of certain 
molecules.^*' In this paper, we will show how a disconti- 
nuity of the XC kernel enters when determining the Fukui 
function and that this quantity is in fact what gives the 
largest contribution in describing the reaction toward a 
nucleophilic reagent. 

The paper is organized as follows. In Sec. II we start 
by deriving an expression for the Fukui function in terms 
of the XC kernel. In Sec. Ill we discuss discontinuities 
in DFT with a particular emphasis on the discontinu- 
ities of the XC kernel and show their importance for 
the formula derived in Sec. II. A numerical investiga- 
tion in terms of two-electron model molecular systems in 
the exact-exchange (EXX) approximation is given in Sec. 
VI. Finally, we give our conclusions in Sec. V. 



II. FUKUI FUNCTION 

In order to define quantities such as the chemical po- 
tential or Fukui functions which involve derivatives with 
respect to the number of particles the theory must be 
generalized to systems which involve fractional charges. 
This implies an ensemble description in terms of states 
with different electron numbers. For an average number 
of electrons N = iVo 4- w, where TVg is an integer and 
< cli < 1 Perdew et al.^^ proposed an ensemble of the 
form 



7^ = (1 -'^)l*Wo)(*A^ol +^i*Wo + l)(*Wo 



-l-ll 



(1) 



where '^k is the ground state wave function of k particles. 
Similarly, for = iVo — 1 + cj we can define 



7< = {l-u)\^No-l){-^No-l\+^\^No){'^No 



(2) 



Using these ensembles any derivative with respect to A^ 
is equal to the derivative with respect to lo. 

The ensemble ground-state energy E{N) will consist of 
straight line segments between the values at the integers. 
Hence, the chemical potential /j, = dE/dN , i.e., the slope, 
on the — /+ side of E{Nq) is equal to the negative of the 
ionization energy and affinity (I/A), respectively. 
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The Fukui function is defined as 



fir) 



dn{r) Sfi 



dN 5w{y) ' 



(3) 



where n(r) is the electron density and w{y) the external 
potential. These two definitions are equal which is easily 
seen from the identity 6E/6w{r) ~ n{r). Using the first 
definition we immediately see that the Fukui function 
is constant (independent of N) between integers. The 
limits N is, however, usually what is of interest 

and we therefore write 



/+(r) = nNo+i{r) - nN„ir) 



(4) 
(5) 



where here and in the following the +/— sign refers to 
these different limits. These results are a direct con- 
sequence of using the above defined ensembles, where 
the species is assumed to be completely independent of 
its environment. To address this aspect, the concept 
of chemical-context dependent reactivity descriptors has 
been proposed. 

So far we have not used DFT but it appears to be a 
natural framework for calculating reactivity indices like 
the Fukui function which is defined in terms of the den- 
sity only. The extension of DFT to fractional charges^^ 
allows us to formally take derivatives with respect to N. 
The density is usually calculated within the KS frame- 
work which assumes the ensemble density to be non- 
interacting ensemble w-representable. The ensemble KS 
potential, 14, has been the subject of several investiga- 
tions all showing that the XC part Wxc(r) must have a 
discontinuous behavior as a function of A^j22,23 topic 
we will elaborate further on in the next section. 

In order to derive expressions for the Fukui function 
in terms of KS quantities we start by writing the density 
corresponding to Eqs. (1-2) in terms of KS orbitals 

No 



>(r) = J2\^tir)\''+c^\VNo+iir)\^ (6) 
fc 

No-l 



We notice that the orbitals depend on lu (or N) since the 
ensemble KS potential will be different for every value of 
Lo. Keeping this in mind we can take the derivative with 
respect to N and determine the Fukui functions^^ 



N„ 



f-ir) = l^.o.i(r)P+E%i^ (8) 



No 



r(r) = \^Noir)\'+J2 



dN 



dN 



(9) 



The superscript uj on the orbitals is now dropped since 
the limit N N^ has been taken. The orbitals are 



continuous with respect to A'' and are therefore unam- 
biguously determined by the A^o-system. The frontier 
molecular orbital (FMO) approximation corresponds to 
ignoring orbital relaxations, i.e., f~^{r) « |(p7Vo+i(r)P 
and /^(r) « \(pis[g{r)\'^ . In the formulation by Parr et 
al.^'^ orbital relaxations due to Coulomb interactions can 
be taken into account via the second terms on the right 
hand side of Eqs. (8-9). The differentiation of the or- 
bitals can easily be performed using the chain rule 



dN 



dr 



6Vs{r') ON ■ 



(10) 



A variation in the number of particles will induce a vari- 
ation in the KS potential via the density 

dVsjr') _ /-^„„ ^bH(rO+t;.c(rO] ,,„„, 
dN ^ J 6n{r") ' 

= J dr"[t-(r',r") + /xc(r',r")]/(r"), (11) 

where WH(r) = / dr''i;(r, r')n(r') is the Hartree potential, 
V is the bare Coulomb interaction, or the Hartree kernel, 
and /xc is the XC kernel. It is then easy to see 



/±(r) = /±(r) 



J dr'rfr"x.(r,r')/±,,(r'. 



')/^(r"),(12) 



where Xs = Sn/dVs is the KS density response function, 
/hxc = w + /xc and Jq is the Fukui function in the FMO 
approximation defined above. Here we have been careful 
in taking the limit A^ — ?■ A^^*^ of the XC kernel since it 
has been recently shown that /xc has discontinuities.^^ 
Defining the matrix 

K^{v,v')^5{v,v')- j dr"x.(r,r")/±xc(r",r') (13) 

we can recast Eq. (12) into 

/±(r) = j dr'[/^±(r,r')]-Vo^(r')- (14) 



This formula was obtained by Cohen et al 



26,27 



apart 



from the fact that we here allow the kernel to have dis- 
continuities. That these play a dominant role when cal- 
culating the Fukui function /+ will be shown in Sec. IV. 

We will now show an alternative derivation of Eq. (14) 
based on the equivalent definition of the Fukui function 
as the functional derivative of the chemical potential with 
respect to the external potential. We begin by eval- 
uating the chemical potential and show its equivalence 
with the highest occupied eigenvalue of the KS system. 
The ground-state energy is a functional of the density 
and can be written as 



E[n] = Ts[n] 



1 



(irdr'ri(r)u(r, r')ri(r') 



j drw{r)n{r) + 



(15) 
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where Ts is the non-interacting kinetic energy functional 
and i?xc is the XC energy. Taking the derivative with 
respect to the number of particles we write 

l|-§+/^^Hr)+.H(r)]/(r) + ^, (16) 
I 



where we have identified the Fukui function. The deriva- 
tive of the kinetic energy can easily be evaluated once 
written in terms of occupied KS eigenvalues e^. Let us 
first focus on an ensemble of the form of Eq. (1), i.e., 
with N = Nq + Lo. We can then write 



ON 



d 

dN 



'No 

E 



Ek + uJENo+i ^ J '^^ [^('") + ^H(r) -i- Wxc(r)] n(r) 



dr [■w{r) + UH(r)] /(r) -|- 



ON 



(17) 



The eigenvalues are functions of N via the KS potential 
and we find after a few manipulations 



= £jVo+i - / drwxc(r)/(r) + 



dN 



where we have used the identity 



dN 



cJi'Wxc(r)/(r) 



;jVo+i, (18) 



(19) 



and the definition of the XC potential Wxc = 5E^c/ 5n. 
The same steps can be performed for the ensemble in 
Eq. (2) and we find similarly 



dE< 
dN 



£ No- 



rn 



Eqs. (18) and (20) thus prove that the highest occupied 
eigenvalue must be equal to the chemical potcntiaP^ and 
should therefore not change with N. Many problems with 
existing functionals are related to a lack of this straight- 
line behavior when extended to fractional charges. '^^■^'^ In 
the limit N 



Nq we have 



dE 
dN 



= e 



dE 
'dN 



'No 



(21) 



'LUMO 



, i.e., the lowest unoccupied 



In this limit e^Vg+i 
KS orbital obtained from the KS potential in the limit 
N -^N+ {V+). In the same way — ^homc i-^-j the 
highest occupied KS orbital obtained from V~ . These 
different limits are important to keep since a constant 
shift (or discontinuity) in Wxc has been shown to occur 
as an integer is crossed. The discontinuity in Uxc is in 
general positive, shifting the KS affinity As — er 
to the true affinity A = e 



+ 

LUMO 



'LUMO 



in order to obey the 
relation in Eq. (18). Before continuing the discussion 
on discontinuities (Sec. Ill) we will determine the Fukui 
function from 



fir) 



6 dE Se 



Na/Na + l 



5ii 

Sw{r) Sw{r) dN Sw{r) 



(22) 



The variation of an KS eigenvalue with respect to the ex- 
ternal potential can straightforwardly be obtained from 



first order perturbation theory. Taking into account that 
varying w will induce a variation of vn and v^c via the 
density we find 

/±(r) = /±(r) 

+ 1 dr'dr"x(r,r")/±,,(r",r')/o±(r').(23) 



Using the relation 

X = Xs + Xs[v + fKc]x 



(24) 



it is easy to see that Eqs. (23) and (14) are equivalent. 
In the following we will prefer the use of Eq. (23) since 
any discontinuity of /xc enters linearly. 



III. DISCONTINUITIES IN DFT 

It has been shown that the ground-state energy 
exhibits derivative discontinuities at integer particle 
number. For non-interacting electrons in a system with 
discrete energy levels this happens only when N crosses 
an integer for which a new orbital with a different eigen- 
value starts to be occupied. In this case, the effect is 
mainly due to the kinetic energy. In an interacting sys- 
tem with discrete energy levels a large part of the deriva- 
tive discontinuity will be contained in E^c and the dis- 
continuity will even show up when N crosses an integer 
lying within a shell with the same eigenvalue. 

A derivative discontinuity in E^c gives rise to a discon- 
tinuity in Wxc in the form of a constant shift Axe- The 
same discontinuity also leads to discontinuities in the XC 
kernel /xc which arc of more complex nature than those 
of Wxc- In this section wc bricfiy review the discontinuities 
of Wxc and /xc- 



A. XC potential 

Any discontinuity of v^c is related to a derivative dis- 
continuity in E-x^c[n[w, N]]. Using the chain rule we can 
write 



dE^ 
dN 



rfrwxc(r) 



dn{r) 
dN 



(25) 
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This identity has to hold true for any value of A'' and 
we can use it to formally write down the value of the 
discontinuous shift Axe at Nq. For N > Nq wc write 



(25) 



(r) 



, (r) + Axe (r) and insert this expression in Eq. 



dE> 
dN 



J dr [z;x-(r)+Axe(r)]/(r). 



(26) 



Now we take the limit N 
ment we find 



^0+ 



and after a rearrange- 



Axe = 



dN 



- / rfri,- (r)/+(r), 



(27) 



where we have used the fact that the Fukui function inte- 
grates to unity. If dE^c/dN is discontinuous at Nq, Axe 
will be finite. As an example we can use the case where 
i?xc is a functional of the KS Green function Gg and 
SE-x^c/^Gs ~ — iExe, where Sxe is called the self-energy. 
Using Eq. (27) we find the celebrated MBPT formula for 
the discontinuity^^ 



Axe = 



J drdrVL(r) [Exe(r, r', e+^j^o) 
-v-Ar')S{r-r')] ^L(r'), 



(28) 



where (^l is the LUMO orbital. In the next section 
wc will use this formula in the EXX approximation for 
which the self-energy corresponds to the Hartrec-Fock 
(HE) self-energy but evaluated with KS orbitals. 

A discontinuous shift in the XC potential implies also 
a shift in the eigenvalues with the same magnitude. We 
can thus write Eq. (21) in the previous section as 



dE 
dN 



'LUMO + ^xc, 



(29) 



that is, the affinity is equal to the KS affinity plus the 
discontinuity, a well-known result. 



B. XC kernel 

In order to determine the discontinuities of the XC ker- 
nel we start by noting that for particle number conserving 
variations of the density /xc is only defined up to the sum 
of two arbitrary functions gxe(r) -|- gxc(r'). This obser- 
vation follows immediately after inspecting the definition 

of /xe 



SEy,c = y"dr'dr/xc(r. 



•')Sn{r')Sji{r). 



(30) 



When we allow for non particle conserving density vari- 
ations the arbitrariness disappears but leaves the func- 
tional discontinuous. In the case of the XC kernel these 
discontinuities will have to take the form 

/x+e(r', r) - f-{r\ r) = gxe(r) + 5xe(r'). (31) 
In order to determine ^xe with the same procedure used 
for the XC potential we study the quantity 



5 5gxe 

5w{r) dN 



(32) 



Varying Eq. (25) with respect to the external potential 
allows us to write this quantity in terms of the XC kernel. 
The kernel can then be written as /^(r', r) = /x^(r', r) + 
(7xe(r', r). Taking the limit N — )• Nq implies g^di^' , r) — > 
ffxe(r) + 5xc(r') and we arrive at 



(irx(ri,r)gxc(r) 



S dE^ 



6w{r) dN 



J drdr'x(ri,r)/x-(r,r')/+(r') ~J rfrt^+eW^ 



El 

Sw{ri) 



(33) 



From this equation (7xe is only determined up to constant. 
This constant is, however, easily fixed by considering the 
the second derivative of E^c with respect to TV 



dr/+(r)5xe(r) 



d'-E, 



dN^ 



J drdr'/+(r)/x-(r,r')/+(r'), (34) 



yielding a condition to be imposed on Eq. (33). The 
function ^xe obtained via Eqs. (33-34) was recently an- 
alyzed in Ref. 25 showing a diverging behavior of the 



form 




(35) 



as r oo. The diverging behavior can be deduced by 
performing a common denominator approximation to Eq. 
(33).^^ That this is indeed a reliable approximation to Eq. 
(33) was shown in Ref. 25. We can now go back to Eq. 
(23) and add the appropriate term to /^(r) arising from 
the discontinuity of /xe- We find 



5 



/+(r) = /+(r) + J «r'x(r,r')[/-e(r',r") + t'(r',r")]/o+(r") + / rfr'x(r, r')5xc(r'), 



(36) 



or 



/+(r)-/+(r)+ / dr'x(r,r')5xc(r') 



(37) 



which defines /^(r) as the Fukui function when the dis- 
continuity is not taken into account. 



IV. EXX APPROXIMATION 

In order to quantify the relative importance of the two 
different contributions to /^(r) (/^ and g^c) we have 
performed a numerical study on model ID molecular sys- 
tems in the EXX approximation. The EXX functional 
is known to contain the derivative discontinuity at even 
integers when defined on densities corresponding to spin- 
compensated ensembles composed of states with different 
number of particles. 

The EXX energy functional is an implicit functional of 
the density given by'^^ 



- J drdr'j{r,r')v{r, r')j{r,r'), 



(38) 



where the density matrix is given by "f{r, r') = 
^^'^'^ (/jfe(r)(/3fc(r') after spin-summations have been per- 
formed. The corresponding EXX potential = SE^/Sn 
can be evaluated as 



dr' 



SE^ (Sn(r') 
Sn{r') 6Vs{r)' 



(39) 



Since Ey, [7] is an explicit functional of the density matrix 
the derivative 6Ex/SVs{r) is most conveniently evaluated 



using the chain-rule 



SE^ 



SVsir) 



dridr2 



SE^ (57(r2,ri) 
57(r2,ri) SVsir) 



(40) 



The potential is only determined up to constant by Eq. 
(39) but this constant may be fixed using Eq. (25). The 
discontinuity in the EXX potentiaP^ is given by Eq. (41) 



drdr'ipi^lr) [— 7(r, r')w(r, r') 
-v-,{r')5{r~r')]Mr'), (41) 



In the case of two electrons the EXX potential takes 
a particularly simple form being equal to minus half the 
Hartree potential Wx(r) = —1/2 J dr'v{r,r')n{r'). The 
EXX kernel is then also easily evaluated resulting in 
/x(r',r) = — l/2v(r',r). The expressions for the poten- 
tial and kernel are evaluated for — > 2^. In the limit 
TV — > 2+ we have to add the discontinuity Ax to and 
two functions to the kernel 

/x+(r',r) - -^Kr',r)+5x(r)+gx(r'). (42) 

We will now use Eq. (33) to determine g^. Using the 
chain-rule 



Sw{r) 



= / dr' 



SVsir') 



SVs{r') dwir) 



(43) 



this equation can be written in terms of Xs only and we 
find 



drxs(ri,r)5x(r) 



S dE^ 



SVs dN 



i J drdv' Xs(ri,r)v(r, r')|(^L(r') 



dv f^(r) 



5Vs{Tl) 



(44) 



To calculate the Fukui function we notice that only the 
quantity xffx is needed. The arbitrary constant can there- 
fore be left undetermined. We also notice that 

/s+(ri) = ./o+(ri) + ^ j dvdv'x{ruv')v{v\v)f+{v). (45) 

In the following section we will solve these two-electron 
equations for model molecular systems to illustrate the 
importance of the discontinuity. 



A. Numerical results for model systems 



Our model systems consist of ID molecules where the 
singular Coulomb interaction has been replaced with 
a soft- Coulomb interaction with softening parameter 1. 
The inter-particle Coulomb interaction is thus model by 
(2^1 — X2)'^ 4- 1 and the external nuclear potentials by 
Zj \Jx\ -|- 1, where Z is the nuclear charge. This model 
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6 8 10 12 14 



a;[a.u.] 

FIG. 1: The Fukui function of a ID Be^+ atom with the 
nucleus located at 10 a.u.. For comparison also the density 
difference An{x) — nN=3(a;) — nN=2(a;) is plotted. The con- 
tribution from the discontinuity is seen to have the largest 
effect. 




8 10 12 14 

x[a.u.] 



FIG. 2: The Fukui function of a ID He^+Be^+ molecule with 
the nuclei located at 8 and 12 a.u.. For comparison also the 
density difference An{x) = nN=3(2;) — nN=2(a^) is plotted. The 
contribution from the discontinuity is seen to have a large 
effect over the He atom. 

has been used extensively in the literature'^^''^* to mimic 
real 3D diatomic molecules with results in qualitative 
agreement, enough also for our purposes. The equations 
derived in the previous sections arc now all evaluated 
consistently using the model potentials. 

The first system we study is a ID Bc^^ {Z = 4) 
atom. In Fig. 3 we show the Fukui function cal- 
culated using different approximations. The black solid 
curve corresponds to using the FMO approximation, i.e., 
/"•"(r) « /o'(r) = |(pL(r)p. The dotted curve is obtained 
by calculating the densities of the N — 2 and = 3 
systems separately and then subtract them. For a func- 
tional with a linear behavior between the integers this 



approach should give the same result as calculating the 
derivative /"*". The derivative is given by the red long 
dashed curve and we see that although the qualitative 
features now agree there are still some discrepancies. The 
fact that An and /"'" do not agree perfectly is related to 
the well known fact that the EXX functional is not linear 
between integers. The green dashed and blue solid thin 
curves shows the different contributions to the derivative 
(Eq. (36)). It is remarkable to see the effect of the dis- 
continuity, which is very large compared to using only the 
fs approximation. The effect is to give a negative con- 
tribution at the nucleus and slow down the asymptotic 
decay. What is perhaps of more interest is the location 
of the peak which is seen to be shifted compared to the 
FMO peak, and hence becomes in good agreement with 
the An result. 

Next, we turn to a molecular system composed of ID 
He^+Be^"'' (Fig. 2). What is particular about this sys- 
tem is that the HOMO and LUMO orbitals are spatially 
well separated. The HOMO is located at the Be site and 
the LUMO at the He site. The density response func- 
tion contains only products of occupied and unoccupied 
orbitals, which decay exponentially with nuclear separa- 
tion. In this case it is thus clear that the correction in 
fs can only be very small, unless /~ becomes very large, 
which is not the case since /~ = l/2v. The blue solid 
thin curve in the figure confirms this fact. We have, how- 
ever, seen that the discontinuity is a diverging function 
and can thus compensate for small excitation functions 
in the response function. Indeed, we find a large con- 
tribution to the Fukui function from the discontinuity. 
Also in this case we see that the peak position is im- 
proved compared to the FMO approximation. From an 
LDA-type of functional such improvement could not be 
achieved due to the lack of a diverging discontinuity. This 
further suggests that in these cases /"*" would be more ac- 
curately calculated by subtracting densities at different 
integer particle numbers. 

To test the consistency of our results we have also ob- 
tained the Fukui function from the numerical derivative, 
i.e., we have calculated 

fNderf \ _ nN=2+AN{x) - nN=2{x) . 

/ {X) - l4Dj 

for small values of AiV using the ensemble of Eq. (2) in 
the EXX functional. When AA^ this result should 
coincide with /+ calculated from Eq. (37). In Fig. 3, we 
show that they, indeed, agree very well. 

V. CONCLUSIONS 

In this paper we have derived equations for the Fukui 
function of DF CRT using a KS reference system. Our 
central result is that a discontinuity of the XC kernel 
enters when calculating the Fukui function for a nucle- 
ophilic attack (/"*"). The importance of this result has 
been demonstrated in model molecular systems, where 
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X [a.u.] X [a.u.] 



FIG. 3: The difference between the analytic {f^) and the 
numerical {n{N = 2 + AN) - n{N = 2)} /AN derivatives at 
AN = 0.1,0.01,0.001,0.0001. To the left: Re'+Be'^+ and to 
the right: Be^+. 

the effect of orbital relaxation has been shown to be al- 
most entirely due to the discontinuity. From this we can 
conclude that in any system where orbital relaxations arc 
important the discontinuity must be incorporated. These 
conclusions are of course based on model systems but 



they naturally carry over to real three dimensional atoms 
and molecules. An implementation of the full EXX func- 
tional for a molecule is, however, quite demanding and 
we therefore leave such investigation for future work. 

Wc would also like to point out that using the HF 
method to calculate Fukui functions'^^ in the FMO ap- 
proximation as opposed to the KS method might yield 
very different results for /+. In HF theory the orbitals 
are determined from a non-local potential, yielding very 
different virtual orbitals. If the HF LUMO orbital is 
closer to the true Fukui function is, however, hard to say. 

As a final remark we have in this paper only considered 
the Fukui function. A whole set of other reactivity de- 
scriptors exist and, in general, discontinuities will show 
up in the derivatives when using a KS reference system. 
In particular, we believe that for the calculation of the lo- 
cal hardness or hardness kernels*""'*^ discontinuities will 
be an essential ingredient. 



^ K. Fukui, T. Yonezawa, and H. Shingu, J. Chem. Phys. 

20, 722 (1952). 
^ K. Fukui, Science 218, 4574 (1982). 

^ R. G. Parr and W. Yang, J Am Chem Soc 106, 4049 
(1984). 

* P. K. Chattaraj, A. Cedillo, and R. G. Parr, J. Chem. 

Phys. 103, 7645 (1995). 
^ P. Itskowitz and M. I. Berkowitz, J. Phys. Chem. A 101, 

5687 (1997). 

^ P. W. Ayers and R. G. Parr, J. Am. Chem. Soc. 122, 2010 
(2000). 

R. G. Parr, W. Yang, P. W. Ayers, and M. Levy, Theor. 

Chem. Acc. 103, 353 (2000). 
^ P. W. Ayers, R. C. Morrison, and R. K. Roy, J. Chem. 

Phys. 116, 8731 (2002). 
® P. Geerlings, F. D. Proft, and W. Langenaeker, Chem. Rev. 

103, 1793 (2003). 
^° P. Geerlings and F. D. Proft, Phys. Chem. Chem. Phys. 

10, 3028 (2008). 
" P. W. Ayers and R. G. Parr, J. Chem. Phys. 129, 054111 

(2008). 

D. C. Ghosh and N. Islam, Int. J. Quantum Chem. Ill, 
11 (2010). 

R. G. Parr and W. Yang, Density- Functional Theory of 
Atoms and Molecules (Oxford University Press, New York, 
1989). 

P. Hohenberg and W. Kohn, Phys. Rev. B864, 136 (1964). 
U. von Earth, Physica Scripta T109, 9 (2004). 
W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 
J. P. Perdew, in Density Functional Methods in Physics, 
edited by R. M. Dreizler and J. da Providencia (Plenum, 
New York, 1985). 

C.-O. Almbladh and U. von Earth, in Density Func- 
tional Methods in Physics, edited by R. M. Dreizler and 
J. da Providencia (Plenum, New York, 1985). 
J. P. Perdew, R. G. Parr, M. Levy, and J. L. Ealduz, Phys. 
Rev. Lett. 49, 1691 (1982). 
^° L. J. Eartolotti and P. W. Ayers, J. Phys. Chem. A 109, 



1146 (2005). 

M. H. Cohen and A. Wasserman, J. Phys. Chem. A 111, 
2229 (2007). 

P. Gori-Giorgi and A. Savin, Int. J. Quantum Chem. 109, 
2410 (2009). 

E. Sagvolden and J. P. Perdew, Phys. Rev. A 77, 012517 
(2008). 

P. W. Ayers, W. Yang, and L. J. Eartolotti, in Chemical 
Reactivity Theory: A Density Functional View, edited by 
P. K. Chattaraj (CRC Press/Taylor & Francis, 2009). 
M. Hellgren and E. K. U. Gross, In Press, 
ArXiv:1108.3100v2 (2012). 

M. H. Cohen, M. V. Gandugiia-Pirovano, and J. Ku- 

drnovsky, J. Chem. Phys. 101, 8988 (1994). 
^'^ M. H. Cohen, M. V. Gandugiia-Pirovano, and J. Ku- 

drnovsky, J. Chem. Phys. 103, 3543 (1995). 

P. W. Ayers, F. D. Proft, A. Eorgoo, and P. Geerhngs, J. 

Chem. Phys. 126, 224107 (2007). 
23 N. Sablon, F. D. Proft, P. W. Ayers, and P. Geerhngs, J. 

Chem. Phys. 126, 224108 (2007). 
^° T. Fievez, N. Sablon, F. D. Proft, P. W. Ayers, and 

P. Geerhngs, J. Chem. Theor. Comp. 4, 1065 (2008). 

J. P. Perdew and M. Levy, Phys. Rev. E 56, 16021 (1997). 
^2 A. Cohen, P. Mori-Sanchez, and W. Yang, Science 321, 

792 (2008). 

P. Mori-Sanchez, A. J. Cohen, and W. Yang, Phys. Rev. 
Lett. 102, 066403 (2009). 

Notice that there are cases where dE^c/dN is discontin- 
uous, but where this discontinuity is exactly canceled by 
the discontinuity of the Fukui function (see e.g. the Hartree 
functional). In this case. Axe is thus zero. 
J. E. Krieger, Y. Li, and G. J. lafrate, Phys. Rev. A 45, 
101 (1992). 

T. Grabo, T. Kreibich, S. Kurth, and E. K. U. Gross, in 
Strong Coulomb correlations in electronic structure calcu- 
lations: beyond the Local Density Approximation, edited 
by V. Anisimov (Gordon and Breach, Amsterdam, 2000). 
D. Tempel, T. J. Martnez, and N. T. Maftra, J. Chem. 



8 



Theor. Comp. 5, 770 (2009). 

N. Helbig, I. V. Tokatly, and A. Rubio, J. Chem. Phys. 
131, 224105 (2009). 

R. Balawender and P. Geerlings, J. Chem. Phys. 123, 
124103 (2005). 

P. W. Ayers and R. G. Parr, J. Chem. Phys. 128, 184108 



(2008). 

M. Berkowitz and R. G. Parr, J. Chem. Phys. 88, 2554 
(1987). 

P. Chattaraj, D. R. Roy, P. Geerlings, and M. T. Sucarrat, 
Theor. Chem. Acc. 118, 923 (2007). 



